// @HEADER
//*********************************************************************//
//  SiYuan: A numerical PDE solver                                     //
//  Copyright (2022) YUAN Xi                                           //
//  This Software is released under the BSD 2-Clause license detailed  //
//  in the file "LICENSE" in the top-level SiYuan directory            //
//*********************************************************************//
// @HEADER

#ifndef _EIGEN_OBSERVER_WRITE_TO_EXODUS_HPP
#define _EIGEN_OBSERVER_WRITE_TO_EXODUS_HPP

#include "Teuchos_RCP.hpp"
#include "Teuchos_Assert.hpp"

#include "Tpetra_Map.hpp"
#include "Tpetra_CrsGraph.hpp"
#include "Tpetra_CrsMatrix.hpp"
#include "Tpetra_Import.hpp"
#include "Tpetra_Export.hpp"

#include "Panzer_STK_Interface.hpp"
#include "Panzer_GlobalIndexer.hpp"
#include "Panzer_STK_ResponseEvaluatorFactory_SolutionWriter.hpp"
#include "Panzer_STK_Utilities.hpp"

#include "AnasaziTypes.hpp"
#include "Thyra_TpetraThyraWrappers.hpp"

namespace SiYuan {

template< typename ScalarT >
class EigenObserver_WriteToExodus {

    typedef Tpetra::MultiVector<ScalarT>      TMV;

public:
    
  EigenObserver_WriteToExodus(const Teuchos::RCP<panzer_stk::STK_Interface>& mesh,
				  const Teuchos::RCP<const panzer::GlobalIndexer>& dof_manager,
				  const Teuchos::RCP<const panzer::LinearObjFactory<panzer::Traits> >& lof,
          const Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > & response_library) :
      m_mesh(mesh),
      m_dof_manager(dof_manager),
      m_lof(lof),
      m_response_library(response_library)
  { 
      // get all element blocks and add them to the list
      std::vector<std::string> eBlocks;
      mesh->getElementBlockNames(eBlocks);

      panzer_stk::RespFactorySolnWriter_Builder builder;
      builder.mesh = mesh;
      m_response_library->addResponse("Main Field Output",eBlocks,builder);
  }
    
  void observeSolution(Anasazi::Eigensolution<ScalarT,TMV>& sol)
  { 
    Teuchos::RCP<TMV> evecs = sol.Evecs;
   // Teuchos::RCP<const Thyra::MultiVectorBase<ScalarT> > thyra_sol = Thyra::createConstMultiVector(evecs);
    for( int i=0; i<sol.numVecs; i++ )  
    {
      Teuchos::RCP<Tpetra::Vector<ScalarT>> a = evecs->getVectorNonConst(i);
      auto solution = Thyra::createVector(a);
    //  Teuchos::RCP<const Thyra::VectorBase<ScalarT> > solution = thyra_sol->col(i);
      // initialize the assembly container
      panzer::AssemblyEngineInArgs ae_inargs;
      ae_inargs.container_ = m_lof->buildLinearObjContainer();
      ae_inargs.ghostedContainer_ = m_lof->buildGhostedLinearObjContainer();
      ae_inargs.alpha = 0.0;
      ae_inargs.beta = 1.0;
      ae_inargs.evaluate_transient_terms = false;

      // initialize the ghosted container
      m_lof->initializeGhostedContainer(panzer::LinearObjContainer::X,*ae_inargs.ghostedContainer_);

      {
         // initialize the x vector
        const Teuchos::RCP<panzer::ThyraObjContainer<ScalarT> > thyraContainer
            = Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<ScalarT> >(ae_inargs.container_,true);
        thyraContainer->set_x_th(solution);
      }

      m_response_library->addResponsesToInArgs<panzer::Traits::Residual>(ae_inargs);
      m_response_library->evaluate<panzer::Traits::Residual>(ae_inargs);
      
      m_mesh->writeToExodus(sol.Evals[i].realpart);
    }
  }
    
protected:

    Teuchos::RCP<panzer_stk::STK_Interface> m_mesh;
    Teuchos::RCP<const panzer::GlobalIndexer> m_dof_manager;
    Teuchos::RCP<const panzer::LinearObjFactory<panzer::Traits> > m_lof;
    Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > m_response_library;
};

}

#endif
